Modeling river delta formation 
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A new model to simulate the time evolution of river delta formation 
process is presented. It is based on the continuity equation for 
water and sediment flow and a phenomenological sedimentation/ 
erosion law. Different delta types are reproduced using different 
parameters and erosion rules. The structures of the calculated 
patterns are analyzed in space and time and compared with real 
data patterns. Furthermore our model is capable to simulate the 
rich dynamics related to the switching of the mouth of the river 
delta. The simulation results are then compared with geological 
records for the Mississippi river. 
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Introduction 

The texture of the landscape and fluvial basins is the product of 
thousands of years of tectonic movement coupled with erosion and 
weathering caused by water flow and climatic processes. To gain 
insight into the time evolution of the topography, a model has to in- 
clude the essential processes responsible for the changes of the land- 
scape. In Geology the formation of river deltas and braided river 
streams has been studied since a long time describing the schematic 
processes for the formation of deltaic distributaries and inter-levee- 
basins (H 121 [21 mis). Experimental investigation of erosion and de- 
position has a long tradition in Geology [6^. Field studies have been 
carried out for the Mississippi Delta |7l|8l|9l[Tol, the Niger Delta 
llllll2|[T3l . or for the Brahmaputra Delta |I14!. Laboratory experi- 
ments have also been set up in the last decades for quantitative mea- 
surements 115,16,17 18, 19, 20 1 . For instance, in the experimental 
EarthScape (XES) project the formation of river deltas is studied on 
laboratory scale and different measurements have been carried out 
(2Tl|22l|^. 

Nevertheless modeling has proved to be very difficult as the sys- 
tem is highly complex and a large range of time scales is involved. To 
simulate geological time scales the computation power is immense 
and classical hydrodynamical models cannot be applied. Typically 
these models are based on a continuous ansatz (e.g., shallow water 
equations) which describes the interaction of the physical laws for 
erosion, deposition and water flow I24II25| |26 ,27 28 , 29 1. The result- 
ing set of partial differential equations are then solved with boundary 
and initial conditions using classical finite element or finite volume 
schemes. Unfortunately none of these continuum models is able to 
simulate realistic land-forms as the computational effort is much too 
high to reproduce the necessary resolution over realistic time scales. 
Therefore in the last years discrete models based on the idea of cel- 
lular automata have been proposed I30II3 II [32l 1331 l34l 1351 . These 
models consider water input on some nodes of the lattice and look 
for the steepest path in the landscape to distribute the flow. The sed- 
iment flow is defined as a nonlinear function of the water flow and 
the erosion and deposition are obtained by the difference of the sed- 
iment inflow and outflow. This process is iterated to obtain the time 
evolution. In contrast to the former models, these models are fast and 
several promising results have been obtained, but as they are only 
based on flow, a well defined water level cannot be obtained with this 



ansatz. 

Here we introduce a new kind of model where the water level and 
the landscape are described on a lattice grid coupled by an erosion 
and sedimentation law. The time evolution of the sediment and water 
flow is governed by conservation equations. The paper is organized 
as follows. After an overview on the different types of deltas and 
their classification the model is introduced and discussed in details. 
The analysis of the model results and a comparison with real land- 
forms are provided. According to different parameter combinations 
different delta types can be reproduced and interesting phenomena in 
the time evolution of a delta such as the switching of the delta lobe 
can be observed. Finally the scaling structure of the delta pattern is 
analyzed and compared with that obtained from satellite images. 

Classification 

The word "delta" comes from the Greek capital letter A and can be 
defined as a coastal sedimentary deposit with both subaerial and sub- 
aqueous parts. It is formed by riverborne sediment which is deposited 
at the edge of a standing water. This is in most cases an ocean but 
can also be a lake. The morphology and sedimentary sequences of 
a delta depend on the discharge regime and on the sediment load of 
the river as well as on the relative magnitudes of tides, waves and 
currents [361. Also the sediment grain size and the water depth at 
the depositional site are important for the shape of the deltaic de- 
position patterns lifl 1361 [37l 1381 . This complex interaction of dif- 
ferent processes and conditions results in a large variety of differ- 
ent patterns according to the local situations. Wright and Coleman 
IT] |36|[39ll401 described depositional facies in deltaic sediments and 
concluded that they result from a large variety of interacting dynamic 
processes (climate, hydrologic characteristics, wave energy, tidal ac- 
tion, etc.) which modify and disperse the sediment transported by 
the river. By comparing sixteen deltas they found that the Missis- 
sippi Delta is dominated by the sediment supply of the river while 
the Senegal Delta or the Sao Francisco River Delta are mainly dom- 
inated by the reworking wave activities. High tides and strong tidal 
currents are the dominant forces at the Fly River Delta. 

Galloway [41 1 introduced a classification scheme where three 
main types of deltas are distinguished according to the dominant 
forces on the formation process: river-, wave- and ftde-dominated 
deltas. This simple classification scheme was later extended I37II38I 
,42] including also grain size and other effects. 

At the river-dominated end of the spectrum, deltas are indented 
and have more distributaries with marshes, bays, or tidal flats in the 
interdistributary regions. They occur when the stream of the river and 
the resulting sediment transport is strong and other effects such as re- 
working by waves or by tides are minor I36II39I . These deltas tend 
to form big delta lobes into the sea which may have little more than 
the distributary channel and its levee is exposed above the sea level. 
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Due to their similarity with a bird's foot, they are often referred in the 
literature as "bird foot delta" like in the case of the Mississippi River 
Delta 1 36 1 . When more of the flood plain between the individual dis- 
tributary channels is exposed above the sea level, the delta displays 
lobate shape. Wave-dominated delta shorelines are more regular, as- 
suming the form of gentle, arcuate protrusions, and beach ridges are 
more common (e.g., like for the Nile Delta or Niger Delta 1121 1431 ). 
Here the breaking waves cause an immediate mixing of fresh and salt 
water. Thus the stream immediately loses its energy and deposits all 
its load along the cost. Tide-dominated deltas occur in locations of 
large tidal ranges or high tidal current speeds. Such a delta often 
looks like a estuarine bay filled with many stretched islands parallel 
to the main tidal flow and perpendicular to the shore line (like e.g., 
the Brahmaputra Delta). Using the classification of Galloway |41 1 
the different delta types of deltas can be arranged in a triangle where 
the extremes are put in the edges (see Fig. [TJ. 

The Model 

The model discretizes the landscape on an rectangular grid where the 
surface elevation Hi and the water level Vi are assigned to the nodes. 
Both Hi and Vi are measured from a common base point, which is 
defined by the sea level. On the bonds between two neighboring 
nodes i and j, a hydraulic conductivity for the water flow from node 
i to node j is defined as 
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As only surface water flow is considered, Oij is set larger than zero 
only if the water level of the source node is larger than the topogra- 
phy, which means that water can only flow out of a node where the 
water level is above the surface. The relation between the flux Uj 
along a bond and the water level is given by 



hj=o.,(y,-Vi). 



[2] 



Furthermore water is routed downhill using the continuity equation 
for each node 

V^ - VI 
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[3] 



where the sum runs over all currents that enter or leave node i and 
Vi' is the new water level. The boundaries of the system are chosen 
as follows: On the sea side the water level on the boundary is set 
equally to zero and water just can flow out of the system domain. On 
the land the water is retained in the system by high walls or choos- 
ing the computational domain for the terrain such that the flow never 
reaches the boundary. Water is injected into the system by defining 
an input current 7o at the entrance node. 

The landscape is initialized with a given ground water table. 
Runoff is produced when the water level exceeds the surface. The 
sediment transport is coupled to the water flow by the rule, that all 
sediment that enters a node has to be distributed to the outflows ac- 
cording to the strength of the corresponding water outflow. Thus the 
sediment outflow currents for node i are determined via 
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where the upper sum runs over all inflowing sediment and the lower 
one over the water outflow currents. A sediment input current so is 
defined in the initial bond. 



The sedimentation and erosion process is modeled by a phe- 
nomenological relation which is based on the flow strength lij and 
the local pressure gradient imposed by the difference in the water lev- 
els in the two nodes Vi and Vj . The sedimentation/erosion rate dSij 
is defined through 

dS,j = Cl(r - \Uj\) + C2iV* - \Vr ~ V,\), [5] 

where the parameters 7* and V^* are erosion thresholds and the co- 
efficients ci resp. C2 determine the strength of the corresponding 
process. The first term c\(I* — \Iij\) describes the dependency on 
the flow strength lij |44| and is widely used in geomorphology, while 
the second term C2{V* — \Vi — V^|) relates sedimentation and ero- 
sion to the flow velocity, which in the model can be described by 
lijjoij \yi — Vj\. The two terms of [s] are not linearly dependant 
on each other as one may think first by looking at Eq. |2] In fact due 
to Eq[T]there is a nonlinear relation between V and / which leads to 
different thresholds in the pressure gradient and the current. 

The sedimentation rate dSij is limited by the sediment supply 
through Jij, thus in the case dSij > Jij the whole sediment is de- 
posited on the ground and Jij is set to zero. In the other cases Jij 
is reduced by the sedimentation rate or increased if we have ero- 
sion. The erosion process is also supply limited which means that 
the erosion rate is not allowed to exceed a certain threshold T; so 
if dSij < T, then dS'ij — T. Note that in the case of erosion dSij 
is negative. Due to erosion or deposition, the landscape is modified 
according to 
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where the sediment deposits equally on both ends of the bond. The 
new topography is marked with 77-. The same formulae (Eqs. |6]|7} 
also hold in the case of erosion when dSij is negative. 

Iterating Eqns. [T] to |7] determines the time evolution of the sys- 
tem. Finally in a real system subaquateous water currents lead to a 
smoothening of the bottom which is modeled by the following ex- 
pression 

' ^ [8] 

N.N. 

where e is a smoothening constant determining the strength of the 
smoothening process. The sum runs over all nearest neighbors of 
node i. 

Simulation 

The simulation is initialized with a valley on a rectangular N x N 
lattice with equal spacing grid as shown in Fig. |2] The valley runs 
downhill with slope S along the diagonal of the lattice and the hill- 
slopes of the valley increase from the bottom of the valley sidewards 
according to a power law with exponent a. In the simulation shown 
in Fig. [3]the value of a was chosen to be 2.0. Under the sea the land- 
scape is flat with a constant slope downhill. Furthermore we assume 
the initial landscape to have a disordered topography by assigning 
uniformly distributed random numbers to 77^. This variable is then 
smoothed out according to Eq. |8] The water level Vi of the system is 
initialized with a given ground water table. In reality the distance of 
the ground water to the surface is minimal on the bottom of the valley 
and increases uphill. This is obtained in the simulation by choosing 
the water level Vi in an incline plane S below the bottom of the valley. 
The slope of the plane is the same as the slope of the valley S. This 
also keeps the river close to the bottom of the valley. As we are only 
interested in studying the pattern formation at the mouth of the river. 
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the braiding conditions of the upper river only determine the width 
of the delta front. On the seaside when Hi < the water level is a 
constant and set to zero. A sketch of the initial landscape is shown in 
Fig.|2] 

An initial channel network is created by running the algorithm 
without sedimentation and erosion until the water flow reaches a 
steady state. Then the sedimentation and erosion procedure is 
switched on and the pattern formation at the mouth of the river is 
studied. 

According to the dominance of the different processes, com- 
pletely different coastline shapes can be observed. The smoothen- 
ing procedure Eq. |8] leads to the formation of an estuary by re- 
working the coastline at the river mouth, while the stream dominant 
erosion term ci(/* — in Eq. [s] favors the formation of river- 

dominated birdfoot shaped delta. In contrast to this, the second term 
C2 ( V* — I Vi — V, i ) in Eq. [s] which depends on the pressure gradient 
represented by the height difference of the water levels in the nodes 
i and j, produces more classical shaped deltas with several islands 
and channels. These patterns are similar to the distributary structure 
of the Lena or the Mahakam river delta, which are more sea or wave- 
dominated. This difference can be explained by the fact that the first 
term sediments along the main current stream, whereas the second 
term distributes the sediment more equally to the neighboring nodes. 

Figures[3|a-c) show some snapshots of the time evolution of the 
simulation of a birdfoot delta (C2 — 0). A map of the Mississippi 
River is given in Fig. |4]for comparison. In both cases one can see 
how the main channel penetrated into the ocean depositing sediment 
mainly on its levee sides. When the strength of the main channel de- 
creases, side channels start to appear breaking through the sidebars as 
can be seen in the snapshots of Figs. |3|b) and |3jc). At the beginning 
of the delta formation process the sediment transport is equally dis- 
tributed among the different channels and leads to a broader growth 
of the delta front along the coast. With time, the side channels are 
gradually abandoned and the sediment is primarily routed through 
the main channel, thus this dominant channel is growing faster than 
the others forming the typical birdfoot shaped deposits. 

Figure [5j a) shows another type of delta where the smoothening 
of the waves reworks the deposits at the river mouth and distributes it 
along the coast. Here the river could built up only a slight protrusion 
in the immediate vicinity of the river mouth. The same happens in ar- 
eas where the wave currents are dominant and lead to the formation 
of wave-dominated deltas like the Sao Francisco River in Brazil or 
the Nile Delta. A map of the Sao Francisco River Delta is also given 
in Fig. [5|b) for comparison. Here the coast line has been straightened 
by the wave activities and consists almost completely of beach ridges 
which have the typical triangular shape inland. This flattened deposit 
can also be found in our simulation results. As there is no evapora- 
tion included in the simulation, small ponds and abandoned channels 
remain in the sedimented zone instead of disappearing with time. 

Finally, if the term C2{V* — \ Vi — Vj\) dominates the sedimen- 
tation/erosion process a half moon shaped delta with many small is- 
lands and channels appear. This delta type shows more activity in 
the channel network than the others. The channels split and come 
together, and when the main channel blocks its way due to sedi- 
mentation, the whole delta lobe switches to another place. This phe- 
nomenon is called delta switching. During the simulation, the switch- 
ing of the delta occurred several times. 

The best studied delta in the world is that of the Mississippi river 
where the switching of the delta lobes was studied in detail H). The 
switching of the Mississippi Delta during the last 4000 years is well 
documented L7.i8. ,45J . The rich dynamics due to the switching phe- 
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nomenon that is observed in the Mississippi can be also identified in 
our simulations. In the literature Ii36 | three types of switching mech- 
anisms are distinguished. The first type referred as switching type I, 
consists of a lobe switching in which the delta propagates in a series 
of distributary channels. After a certain time, the stream abandons the 
entire system close to the head of the delta and forms a new lobe in an 
adjacent region. Very often this lobe occupies an indentation in the 
coastline between previous existing lobes so that with time the sedi- 
ment layers overlap each other One can find this type of delta switch- 
ing in areas where the offshore slope is extremely low and the tidal 
and wave forces are too small for reworking the lobe 1361 1391 |4T1 . 
In many cases the delta lobes merge with each other forming major 
sheet-type sand banks. This phenomena can be nicely observed when 
comparing the two images of the simulation in Figs.|6[a) and|6jb). In 
the Mississippi Delta this also happened several times in the past and 
the different lobes have today different names. A type I shift of the 
Mississippi delta occurred for example, 4600 years B.C. between the 
Sale-Cypremort and the Cocodrie (4600-3500) Lobe or when the St. 
Bernard Lobe switched to which today is called Lafourche at about 
1000 BC. 

At about 3500 BC the Mississippi river switched far upstream 
from the Cocodrie to the Teche stream tailing a completely new 
course for the river itself and its delta. This type of switching is 
referred as type II switching 1 36] and can also be found in the simu- 
lation. When comparing Figs. |6[b) and|6|c) one can see a major shift 
of the channel far upstream in the deltaic plain so the river takes a 
completely different course and forms a new delta beside. 

The type III of delta switching is referred in the literature as alter- 
nate channel extension 136|. In this case not the complete channel but 
the dominance of sediment flux in one or more distributaries is chang- 
ing with time. This can be described as follows: two or more major 
channels split into several distributaries nearly at the same point at 
the head of the delta. Commonly one of the distributaries is dom- 
inant, so it will carry most of the sediment and water discharge at 
any time. As a result, this active channel will rapidly propagate sea- 
ward, while the other channel will shrivel with time. At some point, 
the slope of the main active channel will decrease and the discharge 
will seek one of the shorter distributaries. With the increased sedi- 
ment flux downstream, the new channel will rapidly propagate into 
the sea. This switching process will repeat several times forming a 
deltaic plain characterized by a series of multiple beach ridges. This 
type of switching can be best observed in the simulation of the bird- 
foot delta in Fig. |3|a-c), where the main path of the sediment flow is 
marked by the red arrows. One can see how side channels emerge and 
are abandoned after a certain time. Nevertheless, a major switching 
of the main channel could not be observed in the simulations. The 
average time between two lobe switchings was found to be around 
1000 years for the Mississippi river [8, 36, 45]. 

At this point, we show that the river delta patterns generated from 
our simulations display geometric features that are statistically simi- 
lar to real river delta structures. More precisely, we analyze the self 
similar behavior of the real and simulated river deltas using the box 
counting algorithm |46 |. The box counting dimension is a quite com- 
mon measure in geomorphological pattern analysis and has been used 
by many authors to characterize river basin patterns and coastlines 

Ezlllllllll. 

For the real satellite picture as well as for the simulated river 
delta, we show in Figs. |7]that the variation with the cell size s of the 
number of cells A*' covering the land follows typical power laws over 
more than 3 decades 

iV~s-°, [9] 
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where the exponent D is the fractal dimension. Moreover, the least 
square fit of this scaling function to the data gives exponents which 
are strikingly close to each other, namely D = 1.81 ± 0.01 for the 
real Lena river delta and D = 1.85 ±0.1 for the simulation. 

Conclusion 

A new model for simulating the formation process of river deltas has 
been presented. It is based on simple conservation laws for water and 
sediment on a lattice grid, coupled by a phenomenological sedimen- 
tation/erosion law. Several interesting features of river deltas like the 
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Sao Francisco „ 

Copper 

Fig. 1 . The figure shows the classification scheme after Galloway TP, where wave- tide- and river-dominated deltas are distinguished in the extremes of the triangle. 
By comparing 16 major river deltas Wright & Coelman | 39 concluded that in the extremes the Mississippi is the most river-dominated delta and the Sao Francisco the 
most wave-dominated one. The delta which is mainly dominated by the tides is that of the Fly river in Papua New Guinea. 




Fig. 2. The figure shows a sl^etch of the initial condition for a simulation. A water current Iq is injected at the upper node and the water levels on the sea boundaries 
are kept constant (Vb = 0). The landscape is initialized as an inclined plane with a disordered topography on the top. The water surface (blue) is parallel to the 
horizontal plane. 




Fig. 3. Time evolution of a birdfoot delta (from left to right). Figure (a) shows the delta after 1.2 million time steps, where the main channel worked into the sea 
depositing sediment mainly on its levee sides. After 2.5 million time steps the main channel has split into two distributaries (b), where the smaller one becomes inactive 
after 5 million steps and a new channel breaks through the sidewalls (c). The main directions of the sediment flow are marked with the red arrows. The simulation was 
run on a 279x279 lattice and the parameters for the water flow were Iq = 1.7 x 10^* and Co- = 8.5. For the sedimentation and erosion the constants were set to 
ci =0.1 and C2 = with a sediment input current of so = 0.00025. The erosion threshold /* was set to /* = 4 x 10^® and the maximal erosion rate was set 
to |T| = 5 X lO^''. Smoothening was applied every 2000 time steps with a smoothening factor of e = 1 x 10~*. The initial depth of the water table at the bottom 
of the valley was set to 5 = 0.0025. 
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Fig. 4. For comparison with ttie simulation results of Fig. [sjthe figure shows part of a map of the mouth of the Mississippi river, where the birdfoot shaped delta can 
be seen clearly. The colors indicate channel deposits , sand ridges , swamps and marshes The figure was generated after f36l . 




Fig. 5. In (a) the simulation of a wave-dominated delta Is shown. While the waves are reworking the coast at the mouth of the river to form an estuary, the river 
deposits sediment and forms large beaches. As the simulation does not include evaporation, the ponds and inactive channels In the deposition zone do not disappear 
as in the map of the real river shown In (b). The parameters in the simulation were A"" = 179, Iq = 1.7 x 10""*, sq = 0.0015, c^r = 8.5, ci = 0, C2 = 0.1 
and /* = 1.3 x 10^*. Smoothening was applied every 200 time steps with a smoothening constant e = 0.01. In (b) we show for comparison a map of the Sao 
Francisco river delta in southern Brazil which is the most wave-dominated delta according to the classification of [41]. The colors In the map (Fig. (b)) indicate channel 
deposits , beach ridges , eolian dunes , marsh-mangroves ^, the floodplain ^ and the uplands . The figure was generated after i36i . 




Fig. 6. The figures show the switching of the delta lobe during the simulation. Comparing Figs, (a) and (b) a type I switching can be identified where the main 
part of the delta lobe is abandoned close to the mouth of the river just before the river splits into several distributaries and forms a new lobe beside. Another type of 
delta switching is shown comparing Figs, (b) and (c) with two snapshots from the simulation. Here the channel switches far upstream and takes a new course to the 
coast forming another delta lobe far away. This is referred as a switching of type II. The parameters for the simulation where Iq = 1.7 x 10^*, sq = 5 x 10^^, 
C2 = 0.0005, ci = and /* = 3.3 x 10^*. The simulation was run on a 179 x 179 lattice with smoothening every 2000 time steps and a smoothening constant 
of e = 0.0001. 



6 I www.pnas.org 
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Fig. 7. Figure (a) shows the scaling behavior of the Lena river delta. On the y-axis we show the logarithm of the number of boxes N{s) ot size s which are necessary 
to cover the subaerial surface is plotted versus the logarithm of the inverse box size. The straight line is a power law fit N ~ with exponent D = 1.81. In the 
inset a satellite picture of the Lena delta is shown. In (b) one can see the scaling behavior of the birdfoot delta from the simulation (c.f. Fig. |3] where the slope was 
calculated to be 1 .85. 
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